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ROUNDOFF-ERROR  ACCUMULATION  IN  ITERATIVE  PROCEDURES 
Robert  T.  Gregory  and  A.  H.  Taub 

I.  Introduction.   In  an  unpublished  paper  (see  appendix)  Taub  and  Wilf 
investigate  stopping  procedures  for  linear  iterative  processes.   In  particular, 
they  study  iterative  methods  for  solving  a  system  of  linear  algebraic  equations 
numerically ,  where  each  method  produces  a  sequence  of  vectors  iu  1  with  the  pro- 
perty that  the  solution  vector  is  approximated  by  the  limit  u   of  the  sequence. 

QO 

They  point  out  that  unless  the  stopping  criterion  is  chosen  carefully,  the  number 
of  iterations  performed  may  be  much  larger  than  necessary  due  to  misleading  in- 
formation given  by  the  cancellation  of  roundoff  error. 

The  mathematical  model  chosen  for  their  investigation  is  extremely  simple 
(see  appendix,  section  II ).   It  is  stated  that  the  single  equation 

(1)  u  .  =  bu  , 

n+1     n' 

where 

(2)  1/2  <  b  <  1, 

is  perfectly  general  if  b  is  thought  of  as  the  largest  eigenvalue  of  the  iteration 
operator  b.  In  this  case 

(3)  u  =  0. 

co 

Machine  calculation  of  (l)  actually  produces 

>M  u  .  =  bu  +  5  . 

n+1     n    n 

inhere  6  is  the  roundoff  error.  Under  the  hypotheses  that  the  5  are  uncorrected 

n 
/andom  variables  with  common  probability  density  function  f  (5),  mean  value  .zero, 

i         .       2 

i-ommon  variance  cf    ,  and  uniform  distribution,  it  is  shown  (see  appendix,  sections 

J.I  and  III)  that  stopping  procedure  Pa  (defined  below)  allows  more  iterations, 

in  general,  than  stopping  procedure  Pp   (also  defined  below).  These  extra 

Iterations,  of  course,  are  wasted. 


One  purpose  of  the  present  investigation  is  to  attempt  to  obtain  machine 
-esults  which  verify  the  theoretical  result  of  the  previous  paragraph  and  thus 
validate  the  assumptions  made  about  the  6  .  However,  results  described  below 
.ndicate  that,  for  the  model  chosen,  the  5  are  not  uniformly  distributed  random 
rariables  when  n  is  large  and  the  iterations  approach  convergence. 

II.  The  Machine  Program.  In  order  to  study  the  matematical  model  mentioned 
ibove,  a  program  was  written  for  the  ILLIAC  to  carry  out  the  computation  indicated 
.n  (h-).     It  is  possible,  for  a  given  value  of  b,  to  try  many  cases,  each  one  with 
i  different  starting  value  un  (generated  as  a  paeudo  random  number)  where  b  satis- 
fies (2).  In  each  case,  for  n  =  1,  2,  . ..,  N,  the  quantities  |u  ,  -  u  |  and 

u  ,=u  I  =  |u  .,  I  are  computed. 

n+1    00  '    '  n+1 '        * 

It  is  shown  in  (9)  of  the  appendix  that  the  variance  of  u  ,  -  u  is  less 
v  '  n+1    n 

;han  the  variance  of  u   .  Hence,  noting  (3),  a  stopping  criterion 

(Pa  );  Stop  the  iteration  when  u  ,-,  -  u  is  "pure  noise," 

should  allow  overiteration  as  Qompared  with  a  stopping  criterion 

(P(3  )%     Stop  the  iteration  when  u  .  -  u   =  u  .  is  "pure  noise."  * 

n+1    00    n+1 

Ihe   difference  in  the  number  of  iterations  using  the  two  criteria  is  given  by 

,l4)  of  the  appendix. 

It  was  anticipated  that  the  quantities  lu  n  -  u  I  and  lu  ,  |  would  remain 

1  n+1    n '     '  n+1 ' 

lonotone  decreasing  until  they  became  small  and  roundoff -error  accumulation  became 

ignif leant,  at  which  time  fluctuations  would  occur.  These  fluctuations,  in  each 

ase,  would  be  detected  and  used  as  a  criterion  for  stopping  the  iteration.  The 

■heory  states  that  the  fluctuation  in  lu  .  -  u  I  should  begin  many  iterations 

1  n+i    n' 

ater  than  the  fluctuation  in  lu  ,  | . 

1    n+11 


'  See  equations    (10)  and    (12)  of  the  appendix  for  the  definition  used  for  the 

alue  of  n  at  which  u     ,    -  u     and  u     ,    are    "pure  noise." 

n+1         n  n+1 

-2- 


Because  of  this  fact,  a  third  quantity  was  computed  by  the  program,  namely 


|u  ,  -  u  +  A  .  where  A  is  a  uniformly  distributed  random  variable  generated 
1  n+1    n    n '  *       n  * 

so  as  to  make  the  variance  of  lu  ,  -  u  +  A  I  equal  the  variance  of  lu  ,1,  that 

1  n+1    n    n '  '  n+1 ' ' 

is,  so  as  to  make  the  fluctuation  in  lu  ,  -  u  +  A  I  begin  at  the  same  time  as  the 

'  >    n+1    n    n '   ° 

fluctuation  in  |u    | ,  on  the  average.   In  other  words,  it  was  hoped  that  by  intro- 
ducing some  random  noise  of  the  form  A  one  might  make  up  for  the  cancellation  of 


roundoff  noise  due  to  the  subtraction  performed  in  forming  |u    -  u  |.  Thus,  a 

stopping  criterion  based  on  |u 

average,  than  the  widely-used  criterion  based  on  |u 

III.  Numerical  Results .  The  machine  actually  normalized  the  quantities  men- 


-  u  +  A  I  would  allow  fewer  iterations,  on  the 
n+1    n    n '  s 


-  u 
n+1    n 


tioned  in  the  previous  section  by  dividing  them  by  |un|.  Since  the  results  of  the 
first  few  iterations  were  of  little  interest,  the  printing  portion  of  the  main 
loop  was  by-passed  except  when  any  one  of  the  three  inequalities, 


(5) 


Un+2 

- 

n+1 

uo 

■uo 

>  0 


(6) 

or 
(7) 

was  satisfied 


u 


h+2 


n+1 


u„ 


u  n  -  u 
n+1    n 


u. 


*  °> 


u  _  -  u  .  +  A  . 
n+2    n+1    n+1 


u. 


u  ,.  -  u  +  A 

n+1    n    11 


u. 


>0, 


Table  I  shows  a  typical  set  of  results,  where  b  =  0.9  and  U  =  O.O966  20^+5  7638. 
The  numbers  printed  should  be  multiplied  by  2'^ ,  since  the  ILLIAC  word  contains 
a  sign  and  thirty-nine  binary  digits,  and  the  results  were  printed  as  integers. 
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It  is  seen  that  u  ,  is  monotone  decreasing  until  it  reaches  a  value  such  that 
n+1 

rounded  multiplication  by  b  produces  the  same  value  again.  When  this  occurs,  the 

differences  become  zero,  and  the  pattern  of  behavior  expected  for  |u    -  u  |  does 

not  appear.  More  than  seventy-five  machine  runs  were  made  with  various  values  of 

id  and  \xn,   and  the  results  obtained  in  each  case  agree  with  the  case  shown0 
0 

Figure  1  shows  the  behavior  of  u    graphically.  Here  we  have  results  for 

n+i 

b  =  0.9375  and  u  =  O.O839  0603  8013.  Notice  that  after  330  iterations  the 

-39 
differences  are  constant  and  equal  to  2  > \     At  this  point  it  is  predictable  what 

u   will  be  for  a  given  u  .  Ultimately,  the  constant  value  u  =  8(2  y)   is 

reached.  Since  the  ILLIAC  rounds  by  adding  2    to  the  product  and  then  truncating, 

we  see'  that 

(8)  (0.9375)(8)(2-39)  +  2'k0   =  (if)(2-56)  +  2~k0, 

=  15(2^°)  +  2-k0, 


=  8(2-59), 


and  u  .  =  u  from  this  point  on. 
n+1    n  * 


IV.  An  Analysis  of  the  Difficulty.   It  was  expected  that  when  u    became 
"pure  noise",  we  would  observe  fluctuations  in  its  value  rather  than  a  leveling-off 
at  a  constant  value.   It  is  this  leveling-off  that  explains  our  failure  to  observe 

the  expected  behavior  of  |u  ,  -  u  I . 

1  n+1    n ' 

One  might  ask  whether  or  not  it  was  accidental  that  u  , ,  in  the  example  of 
Figure  1,  reached  the  critical  value  8(2    ),  and  if  so,  could  this  be  avoided. 


-4- 


TABLE  I 

The  numbers  in  the  last  three  columns  should  he  multiplied 

n 

Un+1 
U0 

Un+rUn 

u   -u  +A 
n+1  n  n 

*0 

uo 

190 

1002             .  U3                 113 

191 

90^            103               113 

193 

726               83                  83 

19k 

65^            73               93 

195 

590            63               63 

196 

526            63               73 

197 

478            51               51 

198 

k26                              51               63 

199 

382            4l               kl 

200 

3^2            kl                                     51 

201 

310            31               ^1 

202 

278            31               31 

203 

250            31               kl 

204 

226            21               31 

205 

206            21               31 

206 

186               21                  21 

207 

166               21                  31 

208 

ikG                             21               31 

209 

13^               11                  11 

210 

126               11                  H 

211 

114               11                  21 

212 

102               11                  11 

213 

9k                              11               11 

214 

82               11                  21 

215 

Ik                                    11                  21 

216 

62               11                  21 

217 

50               11                  21 

218 

k2                                      11                   11 

219 

k2                                        0                    0 

220 

k2                                        0                    0 

221 

k2 

0 

11 

,-39 
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n  is  the  number  of  iterations   u  ,  is  the  approximation 

n+  -39 

to  the  solution  and  should  be  multiplied  by  2    . 


-6- 


The  answer  appears  to  be  negative,  because  after  330  iterations 


(9)  \   =  2M2-Jy), 


n 
and  beginning  at  this  point, 


(10)  un+]_  =  un  -  2  » 

so  that  sixteen  iterations  later 

(11)  un  =  8(2"39), 

as  could  be  predicted.  This  behavior  was  observed  in  every  case  tried, 


Let  us  examine  this  analytically.  Set 

39 


{12)  u  =  Z  .  a.2  , 

x  n   i=l   i   ' 

and  consider  the  case 

(13)  b  =  0.9375 

=  1  -  2-\ 
Then 

W  w  .  h  -  2-^)  H  v-*  +  a-*]  truncated 

.  [u  ♦  2-^°  -  f,  «-*-»  1 1    t  „ 

In         i=l   l      J      truncated 
In  order  to  see  what  the  truncated  value  becomes,  we  write  the  binary  digits 
in  the  following  manner: 


al  a2  a3  ah  a5  a6  '"   a37  a38  a39 


ax  q:2  ...  a33  a^  a^ 


10   0   0 


a36  *37  a38  a39 


The  truncation  takes  place  to  the  right  of  the  thirty-ninth  digit,  following  the 
subtraction. 


■7- 


There  are  three  cases  : 


<36 


Case  II     a56  =  1,  a^:  =  a^Q   =  o^  =  0 

Case  III     a,/-  =  1,  at  least  one  of  Q!,7,  <x,n,  or  a,g  does  riot  vanish, 
In  Cases  I  and  II!  the  truncated  value  becomes 


35      ±  h 
(15)  u  A.  =  u  -  E  -  a,2"  "  , 


n+l    n   i=l  i 


In  Case  III  it  becomes 


3'5 

(16)  u  A.  =  u  -  £  .  a.2-1"4  -  2"39. 

x  '  n+l    n   1=1  i 

Thus  (10)  will  be  satisfied  if 

(17a)  u  =  0.000  ...  OlOxxx, 

(17b)  u  =  0.000  ...  011000, 

'  n  ' 

or 

(17c)  u  =  0.000  ...  OOlxxx, 

n  ' 

that  is  to  say,  whenever 

(18)  8(2"59)  <  un  <  2i+(2"39). 

Once  this  situation  occurs,  it  is  certain  that  the  critical  value  8(2    )  will 
be  reached  during  some  subsequent  iteration. 

V.  Conclusion.   It  appears  from  these  results  that  for  the  model  chosen  the 
assumption  that  the  5  are  uncorrelated,  uniformly-distributed  random  variables 
is  invalid.   It  may  be  that  for  a  more  complex  model  (possibly  one  consisting  of 
a  system  of  two  equations)  the  assumptions  made  about  roundoff  errors  are  valid. 
This  will  be  examined  in  a  subsequent  report „ 
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APPENDIX 
STOPPING  PROCEDURES  FOR  LINEAR  ITERATIVE  PROCESSES 
H.  S.  Wilf  and  A.  H.  Taub 

Introduction 

In  recent  years,  rather  considerable  attention  has  been  paid  to  developing 
more  rapidly  convergent  algorithms  for  solving  the  system 

(1)  x  =  Ax  +  b 

numerically,  where  x,b  are  column  vectors,  and  A  is  an  N  x  N  matrix,   These 

methods  all  have  for  their  object  the  construction  of  an  iterative  scheme, 

(2)  u  n  =  Bu  +  C, 

v  '  n+1     n    ' 

where  the  largest  eigenvalue  of  B  is  as  much  less  than  unity  in  magnitude  as 

possible,  and  such  that  the  solution  x,  of  (l)  can  be  obtained  by  inspection 

from  the  limit  u   of  (2 ) . 
00     v  ' 

The  purpose  of  this  note  is  to  point  out  that  unless  the  stopping  criterion 

for  (2)  is  carefully  chosen,  the  gain  obtained  from  reduction  of  the  largest 

eigenvalue  will  be  entirely  thrown  away  because  of  misleading  information 

given  by  the  cancellation  of  roundoff  error. 

Statement  of  the  Problem 

In  (2),  let  us  take  C  =0,  so  that  u  =0.  This  will  not  affect  the 

'         00 

final  results.  Next,  let  us  treat  the  case  N  =  1,  so  that  B  is  a  scalar  b. 
While  seemingly  drastic,  this  is  actually  perfectly  general  if  b  is  thought 
of  as  the  largest  eigenvalue  of  b,  as  can  be  seen  by  diagonalizing  B  in  (2). 
The  iteration  (2)  then  generates  the  numbers 

(3)  u  =  b  u_. 

n      0 
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What  is  actually  calculated,  however,  is 

(k)  u  ,-  =  bu  +  5  , 

v  '  n+1     n    n' 

where  5  is  roundoff  error,  for  which  we  assume 
n 

HI:  The  6  have  a  common  distribution  function  f (5)  with  finite  variance 

2 
a-     . 

H2:   The  5  are  uncorrelated  random  variables, 
n 

Consider  the  following  two  stopping  procedures: 

(P3  ):  Stop  the  iteration  when  u  ,  -  u  is  "pure  noise." 

(IB  ):   Stop  the  iteration  when  u  -  u  =  u  is  "pure  noise." 
v  K  '  *  n    co   n 

We  see  that  Pa  is  the  one  normally  used  because  the  limit  u   is  un- 

co 

known.  We  will  now  show  that  Pa  will  always  allow  the  calculation  to  go  on 
much  longer  than  FB  would,  the  extra  iterations,  of  course,  being  quite  wasted 
since  the  solution  is  not  being  improved.   This  phenomenon,  as  will  be  seen 
below,  is  directly  the  result  of  "cancellation"  of  noise  in  u  .  -  u  by  the 
subtraction  involved,  which  gives  the  iterates  a  deceptively  significant 
appearance . 
[.  Analysis  of  Pa  and  Pg 

The  solution  of  (k)   is 
(5) 
Thus, 


u  =  bV  +  b11"^  +  bn""2o,  +  ...  +  b8  0  +  5  .  |  . 
n     0        0       1  n-2    n-lj 


(6)    un+1  -  un  =  bn(b-l)u0  +  [bn"1(b-l)50  +  ...  +  b(b-l)6n_2  *  (b-l)5n_1  +  Bj  , 

and  the  variance  is  (since  the  5  are  uncorrelated) 

x  n 


-10- 


(7)  ^(o)  =  (*-i)V2  /jdj^  Lr2, 

o-2      C       ^2n    .   n.2n+l/- 


For  process  Pf3  ,  we  read  directly  from   (5)  the  noise  level 


r%)  =cr2  ji%n 

'l-b 


We  note  immediately  that  for  large  n,  l/2  <  b  <  1, 


(9)  o-l(a)<a{bh 


so  the  noise  level  has  been  reduced  by  the  subtraction. 
Now,  Pec  will  stop  about  when 

(10)  bnuQ  =<rn(cr), 

J2-b2n+b2n+1  '  V2 
1        1  +  b 

Neglecting  b   compared  to  2,   and  solving  for  n  we  find 


2cr2 

|  logV  (l+b)u^ 

(ll)  n    =  ; 

a        log  b 


Process  P£ 

will  stop  when 

(12) 

bn&  =<rn0) 

-fef' 

1/2 
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